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ABSTEIACT 

Smoothed  particle  hydrodynamics  (SPH)  is  a  non-Monte  Carlo  free- 
Lagrangian  method  of  computational  fluid  dynamics  that  has  recently 
expanded  in  scope  beyond  its  parent  application  of  astrophysical 
simulations.  Modelling  terrestrial  flows  involving  boundaries  is  an 
important  extension  of  the  technique.  This  thesis  describes  the 
application  of  SPH  to  model  both  compressible  and  incompressible  flows 
involving  boundaries . 

Simulations  of  shock  wave  reflection  in  a  one -dimensional  shock 
tube  have  been  made.  Different  boundary  methods  have  been  compared  with 
regards  to  their  ability  to  model  a  one -dimensional  reflecting  shock 
wave.  The  simulations  show  that  the  introduction  of  imaginary  particles 
appears  to  be  the  most  effective  method  of  boundary  simulation  for 
transient  shock  reflection.  Simulations  also  show  that  the  choice  of 
artificial  diffusion  treatment  for  best  reflection  is  ambiguous. 

Incompressible  flow  has  been  modeled  using  several  test  problems. 
To  model  incompressible  flow,  modifications  to  the  standard  SPH 
equations  are  necessary.  An  artificial  compressibility  technique  is 
introduced  through  a  fluid  specific  equation  of  state.  Real  viscosity 
has  been  introduced.  Results  are  shown  for  three  test  cases:  Couette 
flow,  Stoke's  first  problem,  and  a  shear  driven  cavity. 
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NOMENCLATURE 

Local  speed  of  sound 
Arbitrary  function 
Diffusive  flux 
Anti -diffusive  flux 
Smoothing  length 

Subscript  defining  individual  particle 
Subscript  defining  a  summation  index 
Mass 

Number  of  particles 
Pressure 

General  variable 
Boundary  position 
Position 

Distance  between  particles 
Time 

Internal  energy 
Velocity 
Weighting  kernel 

Water  equation  of  state  constant 

Artificial  viscosity  constants 
Dirac  delta  function 
Artificial  viscosity 
Viscous  stress 
Density 

Initial  density 
Distance/smoothing  length 
Relaxed  force  distance 
Force  constant 

Diffusion  coefficient  for  FCT 
Monaghan's  viscosity  coefficient 
Water  equation  of  state  exponent 
Antidiffusion  coefficient  for  FCT 
Kinematic  viscosity 
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Chapter  1.  INTRODUCTION 

Smoothed  particle  hydrodjmamics  (SPH)  is  a  free-Lagrangian 
computational  method  introduced  by  Lucy  [1],  and  Gingold  and  Monaghan 
[2].  This  method  was  originally  developed  for  the  simulation  of  such 
astrophysical  phenomena  as  binary  star  interactions  [3],  super  nova 
explosions  [4],  and  the  origin  of  the  moon  [5-7].  Recently,  more 
interest  has  centered  on  the  terrestrial  applications  of  SPH.  Libersky 
and  Petschek  [8]  have  adapted  SPH  to  simulate  high  speed  impacts  of 
solid  objects.  Several  authors  [9-11]  have  simulated  unsteady  gas- 
dynamic  shock  tube  problems.  Henneken  and  Icke  [12]  have  reproduced 
shock  waves  formed  by  steady  supersonic  flow  impinging  on  a  forward 
facing  step.  Other  applications  include  the  modeling  of  volcanic  ash 
flows  [13],  and  simple  free  surface  incompressible  flows  [14]. 

The  basic  premise  of  SPH  is  to  model  a  continuous  fluid  as  a 
discrete  collection  of  interacting  particles.  Each  particle  has 
characteristic  properties  of  mass,  velocity,  position,  and  energy 
associated  with  it.  These  values  are  updated  through  a  timestepping 
algorithm.  Smoothed  particle  hydrodynamics  does  not  need  a  mesh  or 
finite  differences  to  compute  spatial  derivatives.  Derivatives  are 
calculated  through  the  use  of  a  weighting  kernel,  with  particles  loosely 
viewed  as  interpolation  points .  The  kernel  is  a  smoothing  function  with 
a  smoothing  length  set  such  that  the  effects  of  a  central  particle  just 
reach  its  nearest  neighbors.  By  use  of  a  kernel,  the  mesh-tangling 
problems  of  grid-based  Lagrangian  codes  are  avoided. 
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Smoothed  particle  methods  differ  from  Monte-Carlo  methods  because 
they  do  not  use  random  number  generation  to  determine  the  outcome  of 
particle  collisions  or  diffusion.  A  benefit  that  SPH  has  in  common  with 
other  explicit  Lagrangian  methods  is  that  it  is  naturally  suited  for 
massively  parallel  computation.  An  individual  processor  can  be  assigned 
to  an  individual  particle.  As  massively  parallel  machines  become  larger 
and  more  commonly  used,  this  advantage  should  promote  wider  use  of  SPH. 

The  primary  objective  of  this  thesis  is  to  expand  and  evaluate  SPH 
for  situations  involving  solid  boundaries.  Because  astrophysical 
phenomena  typically  have  no  solid  boundaries,  this  area  has  received 
little  attention. 

The  next  chapter  of  this  thesis  introduces  the  SPH  technique.  The 
third  chapter  reviews  work  on  modeling  compressible  flows  involving 
boundaries.  Reflected  shock  waves  in  a  one -dimensional  shock  tube  have 
been  simulated.  Several  different  schemes  have  been  considered  for 
simulating  the  presence  of  solid  reflective  boundaries.  The  effects  of 
an  alternative  artificial  diffusion  treatment  have  also  been  examined. 
The  fourth  chapter  reviews  work  in  extending  SPH  to  model  incompressible 
flow.  This  requires  the  use  of  an  artificial  compressibility  technique. 
Real  viscosity  has  been  added  to  the  equations  of  motion  and  a  driven 
cavity  has  been  simulated.  The  final  chapter  concludes  the  thesis  and 
make  recommendations  for  further  study. 
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Chapter  2.  FUNDAMENTALS 


2.1  Theory 

Smoothed  particle  hydrodynamics  treats  a  fluid  as  a  collection  of 
particles;  each  having  its  own  characteristic  properties  of  velocity, 
energy,  and  mass.  Any  function,  F(r),  can  then  be  estimated  using  an 
integral  interpolant,  with  the  particles  acting  as  interpolation  points: 

<F(r)>  =  jF(r' )W(r-r' ,h)dr'  (1) 

with  W(r-r',h)  being  the  interpolation  kernel  (hereafter  abbreviated  as 
Wij).  If  the  particles  are  distributed  in  the  domain  with  a  number 
density  n(r),  then  Equation  (1)  can  be  approximated  by 

<F(r)>  =  S  Fj/nj  *  (2) 

Here  the  summation  is  taken  over  all  particles.  The  number  density  can 
then  be  represented  by  use  of  the  particles'  characteristic  mass  and 
density  by 

n(rj)  =  p(rj)/mj.  (3) 

Thus  Equation  (2)  can  be  written  as 

<F(r)>  =  2  FjEOj/pj  *  Wi;j.  (4) 

The  equations  of  motion  may  be  derived  from  this  basic  relation. 


4 


2.2  Equations  of  Motion 

2.2.1  Mass  conservation 

Mass  conservation  is  satisfied  by  direct  substitution  of  p  for  F 
in  Equation  (4): 

pi  =  Epjinj/pj  *  Wij  (5) 

or 

Pi  =  SnijWij. 

An  alternative  approach  for  determining  p  is  through  use  of  the  Eulerian 


continuity  equation: 

dp/dt.  +  pV*v  =  0.  (7) 

The  product  of  density  and  the  velocity  divergence  is 

(pV«v)i  =  Vi»(piVi)-ViViPi.  (8) 

Applying  Equation  (4)  yields 

(pV»v)i  =  2(pv)jmj/pj  ViWij  -  SViPjmj/pj  ViWij.  (9) 

Reducing  gives  the  following  form  for  the  divergence  term 

(pV»v)i  =  2(Vj-Vi)mjViWij .  (10) 

The  density  can  now  be  calculated  from: 

dpi/dt  =  -Smj(vi-vj) ‘ViVij.  (11) 


With  this  method  the  density  will  only  change  when  particles  move 
relative  to  one  another.  The  advantage  of  using  the  divergence  form 
(Eq.  11)  in  preference  to  the  mass  summation  form  (Eq.  6)  manifests  as  a 
savings  in  CPU  time.  All  the  subsequent  governing  equations  can  be 
solved  using  only  the  divergence  of  the  kernel ;  the  kernel  itself  needs 
not  be  calculated.  A  disadvantage  of  the  divergence  form  is  that  mass 
is  not  strictly  conserved  [15]. 


2.2.2  Momentum  equation 

The  momentum  equation  for  an  inviscid  fluid  without  body  forces  is 
dv/dt  =  -1/p  ^  VP.  (12) 

It  is  advantageous  to  write  this  equation  in  its  symmetric  form,  so  that 
linear  and  angular  momentum  are  conserved  exactly  [15].  Thus,  we  write 
dv/dt  =  -V(P/p)  -  (P/p^)Vp.  (13) 

Applying  Equation  (4)  yields: 

dVi/dt  =  -  [2(Pj/pj)mj/pj  VWij  +  S(Pi/pi^)pjnij/pj  VWij],  (14) 
Reducing  and  combining  like  terms  gives  the  final  form  of  the  momentum 
equation. 

dvi/dt=  -i:mj*(Pi/p/  +  Pj/p/)ViWij.  (15) 

2.2.3  Internal  energy  equation 

The  rate  of  change  of  internal  energy  is 

du/dt  =  -(P/p)7*v.  (16) 

In  s3TTimetric  form  this  equation  is: 

du/dt  =  -7»(Pv/p)  +  vV(P/p).  (17) 

Applying  Equation  (4)  to  the  energy  equation  yields: 


dui/dt  -  Smj*(Pj/p3*)(Vi-Vj)»ViWij. 


(18) 
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2 . 3  Kernels 

The  kernel  needs  to  satisfy  two  restrictions; 

Xw(r-r' ,h)dr'  =  1  (19) 

and 

limh->o  W(r-r',h)  =  5(r-r').  (20) 

The  kernel  should  also  have  continuous  first  and  second  derivatives. 

One  obvious  kernel  is  a  gaussian  curve ;  another  is  an  exponential 

W  =  I/tt^''"  exp(-(r-r' )Vh^)  (21) 

where  h  is  the  smoothing  length.  This  curve  has  the  appearance  of  a 
gaussian  and  satisfies  the  above  requirements  (Eqs.  19  and  20).  One 
problem  with  this  form  is  that  one  particle  interacts  with  all  of  the 
other  particles  in  the  domain  (albeit  extremely  weakly  with  most).  This 
causes  CPU  time  requirements  to  scale  with  the  square  of  the  nximber  of 
particles . 

A  more  economical  choice  of  a  kernel  is  a  weighting  system  that 
offers  compact  support.  Compact  support  limits  one  particle's 
interactions  to  a  definite  surrounding  volume.  This  allows  the  code  to 
be  written  with  linked  lists  by  superimposing  a  grid  over  the  particles 
which  serves  to  locate  neighboring  particles.  The  grid  needs  not 
conform  to  the  precise  problem  geometry,  nor  is  any  information  stored 
at  nodal  locations;  thus  no  grid  entanglements  are  produced.  The  width 
of  the  grid  cells  is  set  such  that  a  particle  in  one  cell  can  only 
interact  with  particles  in  the  immediately  neighboring  cells.  This 
technique  [17,18]  allows  CPU  time  to  scale  proportionally  with  the 
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number  of  particles.  The  compact  kernel  [19]  used  for  this  paper  is 
defined  by: 

{( 1  - 1 . 5  7'+ .  7  5  7" )  /  ( TTh )  0<7<1 

(.25*(2-7)^)/(^h)  1<7<2  (22) 

0  7>2 

with  7  =  distance/smoothing  length,  ^  is  the  number  of  dimensions  and  o 
is  the  normalization  constant.  In  one,  two,  and  three  dimensions  this 
constant  is  set  to  2/3,  10/77r,  and  I/tt  respectively.  The  smoothing 
length,  h,  is  generally  taken  to  be  twice  the  distance  between 
particles , 

Another  problem  with  the  exponential  form  is  that  information  can 
propagate  over  the  entire  computational  domain  in  a  single  time  step. 

For  application  of  SPH  to  problems  involving  supersonic  flows,  the 
kernel  can  introduce  errors  because  downstream  values  should  have  no 
effect  on  the  flow  upstream  of  the  shocks.  In  SPH  the  kernel  allows 
particles  both  downstream  and  upstream  to  be  used  in  calculating  shock 
values .  A  compact  support  formulation  of  the  kernel  can  help  reduce  the 
downstream  influence. 

2.4  Artificial  Diffusion 

Artificial  diffusion  is  required  to  dampen  nonphysical,  high- 
frequency  oscillations  that  are  inherent  to  Lagrangian  simulations . 
Previously,  artificial  diffusion  has  been  implemented  in  SPH  through 
various  artificial  viscosity  formulations. 
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2.4.1  Monaghan- Gingold  artificial  viscosity 

Monaghan  and  Gingold  [9]  have  formed  the  most  widely  used  form  for 
artificial  viscosity  (hereafter  described  as  M-G  viscosity)  which  is 
implemented  by  rewriting  the  momentvim  equation  as: 

dVi/dt  =  -2mj*(Pi/pi^+Pj/Pj^+nij)ViWij.-  (23) 


where 


Here 


n 


ij 


-  (acij^ij+^//ij^)/Pij  if  Vij*rij>0 

0  if  Vij*ri3<0. 


(24) 


Mir  (hVijTij)/(rij2+r?") ,  (25) 

and  Cij  is  the  local  speed  of  sound.  Artificial  viscosity  is  introduced 
only  when  particles  are  approaching.  The  constants  and  are 

generally  set  equal  to  1,  2,  and  .Olh^  respectively. 


2.4.2  Flux  corrected  transport 

An  alternative  implementation  of  artificial  diffusion,  previously 
unused  in  SPH  simulations,  is  flux  corrected  transport  [20]  (FCT).  Flux 
corrected  transport  is  a  two-step  process  in  which  diffusion  is  applied 
everywhere  in  the  first  (predictor)  stage  to  eliminate  high  frequency 
oscillations,  and  anti -diffusion  is  selectively  added  near 
discontinuities  in  the  second  (corrector)  stage  to  reduce  smearing.  The 
criteria  for  elimination  is  determined  by  the  relative  location  of  the 
particle  with  respect  to  steep  gradients  of  density.  The  diffusion  is 
applied  on  an  individual  particle  basis  and  whenever  a  fluid  quantity  is 
added  to  at  one  point,  it  is  subtracted  from  another.  Therefore,  there 
is  no  net  gain  or  loss  incurred  on  the  system  as  a  whole.  The  general 
algorithm,  following  Book  et  al.  [21],  is: 


1.  Generate  diffusive  fluxes: 

f^j+l/2  =  i^j+1/2  (  ^"j+1 ) 

2.  Generate  anti -diffusive  fluxes: 

~  Mj+i/2(Q  j+i'^  j) 

3.  Diffuse  solution: 

^  j  “  j-1/2 

4.  Calculate  first  differences  of 

=  <?j+i  ■<5j 

5.  Limit  antidiffusive  fluxes; 

S  =  sign  of  (fj+i/2“'*) 

f  1+1/2*^  =  S  max[0,min{SAqj.i/2***,  |fj+i/2*^|,  SAqj+j/j*** )  ] 

6.  Antidiffuse  solution: 

rt  n+1  —  rt  ***  ad  t  fr  ad 

qj  ~  qj  -tj+i/a  +r  j-1/2 

The  diffusion,  u,  and  antidiffusion,  /z,  coefficients  are  dependent  on 
the  computational  scheme.  For  SPH,  u  and  /z  were  optimized  to  be  0.5  and 
0.25  respectively. 

2.5  Time  Stepping 

A  second  order  Runge-Kutta  scheme  was  used  to  march  the  transport 
equations  in  time.  For  simulations  involving  FCT  the  diffusion  and 
anti -diffusion  were  calculated  before  and  after  the  timestep  (no 
calculations  were  made  at  the  intermediate  step) .  Because  SPH  is  an 
explicit  scheme,  there  is  a  stability  constraint  placed  on  the  size  of 
the  time  step.  For  the  shock  tube  problem,  Monaghan  and  Gingold  [9] 
suggested  a  timestep  limit  of  At  <  0.3h/Ci,  where  Ci  is  the  local  speed 
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of  sound  for  the  particle.  All  simulations  were  run  with  a  time  step 
under  this  limit. 

2.6  Boundary  Conditions 

The  treatment  of  rigid  boundaries  has  received  relatively  little 
attention  in  flow  simulations  using  SPH.  There  are  several  potential 
schemes  that  might  be  considered  for  implementing  a  reflecting  boundary: 
imaginary  particles,  repulsive  forces,  immobile  particles,  or  boundary 
modifications  to  the  equations  of  motion. 

Libersky  and  Petschek  [8]  used  imaginary  particles  to  model  the 
impact  of  an  iron  rod  with  a  rigid  surface.  The  imaginary  particle 
scheme  works  by  placing  extra  particles  along  the  outside  edge  of  the 
boundary.  For  every  "real"  particle  within  a  distance  2h  of  the 
boundary,  an  extra  particle  is  placed  symmetrically  on  the  opposite  side 
of  the  boundary.  These  particles  have  the  same  density  and  pressure, 
and  a  velocity  opposite  that  of  their  corresponding  real  particle 
counterpart . 

Because  the  imaginary  particles  must  be  included  in  the 
computation,  this  technique  requires  extra  CPU  time  and  memory. 

Extension  to  multi -dimensional  simulations  may  create  particle  tracking 
problems.  Instead  of  carrying  the  extra  particles  through  the 
summation,  it  is  possible  to  implement  the  boundaries  using  a  series  of 
logical  if-then  statements.  When  a  real  particle  come  within  a 
distance,  2h,  of  the  boundary,  a  series  of  if-then  steps  can  be  used  to 
determine  if  additional  calculations  need  to  be  made.  This  eliminates 
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the  need  to  keep  track  of  imaginary  particles,  allowing  for  easier 
multi -dimensional  coding,  and  reduced  memory  and  CPU  requirements. 

Monaghan  [14]  has  used  the  repulsive  force  technique  for  simple 
free-surface  incompressible  flow  simulations.  Nagasawa  and  Kuwahara 
[13]  have  used  this  technique  for  modeling  volcanic  ash  flows  as  well. 
The  repulsive  force  technique  works  by  introducing  an  interparticle 
force  between  boundary  particles  and  the  approaching  fluid  particles. 

If  a  particle  moves  within  a  preset  distance,  a  repulsive  force  is  added 
to  the  momentum  summation.  This  force  can  be  modeled  as  a  Lennard- Jones 
type  intermolecular  force  [13,14]: 


f(r) 


if  X<Rb 


(26) 


0  if  x>Rb 

where  x  is  usually  set  equal  to  the  initial  particle  separation,  Rt  is 
the  distance  to  the  boundary,  and  ki,  is  determined  by  the  initial 
conditions.  Alternately,  the  repulsive  force  can  be  modeled  by  a  simple 
spring  force : 


f(r)  = 


1 


(27) 


^  K,(x-Rb)*  if  x<Rb 

0  if  x>Rb 

where  x  is  the  relaxed  spring  distance  (taken  to  be  the  initial  particle 
separation) ,  and  <c,  is  the  spring  force  constant  which  is  determined  by 
the  initial  conditions .  Since  this  method  does  not  involve  extra 
particles,  except  for  those  used  to  outline  the  boundary,  there  are  only 
modest  memory  requirements. 

The  immobile  particle,  or  infinite  inertia,  boundary  treatment 


involves  placing  stationary  particles  along  the  boundaries.  These 
particles  are  included  in  the  summations,  and  after  each  iteration  their 
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velocities  are  reset  to  zero.  This  technique  seems  most  applicable  for 
inactive  boundaries  in  which  there  are  no  new  forces  being  applied  to 
the  boundary.  Campbell  [16]  showed  that  in  certain  simulations  where 
the  boundary  pressure  is  known,  boundaries  can  be  implemented  by 
modification  of  the  equations  of  motion.  This  technique  works  for 
problems  involving  a  fluid  being  driven  by  an  external  force  (such  as 
piston  type  or  bullet/gun  barrel  simulations). 
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Chapter  3.  COMPRESSIBLE  FLOW 


3.1  The  Test  Problem 

One  of  the  vehicles  for  testing  the  boundary  and  artificial 
diffusion  treatments  is  a  close -ended,  one -dimens ional ,  adiabatic  shock 
tube.  This  problem  involves  simulating  the  expansion  and  shock  waves 
generated  when  a  wall  separating  high  and  low  pressure  regions  of  a 
perfect  gas  is  suddenly  removed.  Unlike  previous  SPH  work  with  the 
shock  tube  problem,  simulations  were  continued  until  after  the  shock 
wave  reflected  off  the  tube  end-walls.  Inviscid,  one -dimensional  theory 
provides  an  exact  solution  that  is  valid  until  the  reflected  wave 
impinges  upon  the  following  contact  surface.  The  model  shock  tube  is 

1.2  meters  long  and  has  an  initial  pressure  ratio  of  four  to  one.  The 
SPH  simulation  used  450  particles  with  a  width  set  equal  to  the  twice 
the  initial  separation  on  the  low  density  side.  All  simulations  used 
the  imaginary  particle  technique  for  the  boundary  on  the  rarefaction  end 
of  the  tube . 

3 . 2  Discussion 

3.2.1  Pre-reflection 

Figure  1  shows  plots  of  density,  pressure,  temperature,  and 
velocity  in  the  model  shock  tube  0.9  milliseconds  into  the  simulation, 
but  before  reflection.  This  simulation  used  the  M-G  form  of  artificial 
viscosity.  The  SPH  simulation  tracks  the  location  and  magnitude  of  the 
expansion  and  compression  fronts  well.  Notable  errors  include  the  blip 
in  the  pressure  plot,  an  overprediction  in  the  velocity  rarefaction 


Figure  1:  Plots  of  a)  velocity,  b)  pressure,  c)  density,  and  d)  temperature 

0.9  milliseconds  into  the  simulation,  but  before  reflection.  The 
simulation  uses  Monaghan -Gingold  artificial  viscosity  (dashed 
line),  and  is  compared  to  the  exact  solution  (solid  line). 
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wave,  and  smoothing  of  the  shock  wave  fronts.  Monaghan  and  Gingold  [9] 
attribute  the  pressure  blip  to  the  influence  of  the  contact 
discontinuity.  The  smoothing  of  the  shock  fronts  (over  a  distance  3h) 
is  due  to  the  application  of  artificial  diffusion.  These  results  show 
that  our  simulation  matches  previous  SPH  work  [9-11].  The  simulation 
underpredicts  the  temperature  immediately  after  the  shock.  Comparisons 
against  other  temperature  predictions  in  the  SPH  literature  cannot  be 
made  because  these  others  were  either  run  isothermally  [11]  or  did  not 
provide  temperature  results  [9,10]. 

3.2.2  Artificial  diffusion 

The  purpose  of  adding  artificial  diffusion  is  to  damp  out  high 
frequency  oscillations  that  form  in  the  presence  of  flow 
discontinuities,  like  shock  waves.  While  this  is  common  practice  in 
Eulerian  simulations,  the  Lagrangian  nature  of  SPH  makes  it  even  more 
necessary.  Previous  SPH  shock  tube  work  [9,11]  has  tested  a  variety  of 
artificial  viscosity  forms,  and  concluded  that  the  Monaghan- Gingold  [9] 
(M-G)  viscosity  (used  to  produce  the  Figure  1  results)  leads  to  the  most 
accurate  results  for  the  adiabatic  shock  tube  problem. 

Figure  2  shows  plots  of  the  pre- impact  shock  waves  using  flux 
corrected  transport  (FCT) .  Comparison  with  Figure  1  shows  that  both 
artificial  diffusion  techniques  prevent  high  frequency  oscillations  from 
forming.  Flux  corrected  transport  reduces  the  errors  seen  in  the 
rarefaction  overprediction  and  also  reduces  smoothing  of  the  shock 


fronts . 
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Figure  3  shows  the  enlarged  velocity  profiles  predicted  using  M-G 
viscosity  and  FCT,  where  FCT  more  accurately  reproduces  the  shock  front. 
Another  criteria  for  comparison  of  the  two  artificial  diffusion  schemes 
is  the  CPU  time  required  for  calculations.  An  artificial  viscosity 
treatment  requires  a  summation  over  nearest  neighbors,  while  FCT 
requires  only  a  loop  over  the  particles.  Code  using  FCT  ran 
approximately  15%  faster  than  identical  code  using  M-G  viscosity. 

3.2.3  Reflected  shock 

Four  different  end  wall  techniques  were  considered  for  the  shock 
tube  end-walls.  Of  these,  two  methods  either  failed  or  were 
inappropriate  for  the  shock  tube  simulation.  Campbell's  [16]  method  of 
introducing  boundary  conditions  into  the  equations  of  motion  was  not 
appropriate  for  the  shock  tube  problem  because  the  wall  pressure  is  not 
known.  Simulations  using  immobile  particles  failed  because  the  momentum 
of  the  computational  particles  was  sufficient  to  enable  them  to  cross 
over  the  tube  end  walls. 

Figures  4  and  5  shows  density,  velocity,  temperature  and  pressure 
after  reflection  predicted  using  the  imaginary  particle  method.  The 
reflective  boundary  was  modeled  by  placing  fifteen  extra  particles  on 
the  opposite  side  of  the  tube  end.  This  amount  was  determined  by  the 
number  of  real  particles  that  came  within  2h  of  the  boundary  during  the 
simulation.  Figure  4  was  obtained  using  M-G  viscosity,  and  Figure  5 
using  FCT  artificial  diffusion  methods.  Flux  corrected  transport 
produced  nonphysical  oscillations  behind  the  reflected  shock  wave. 

These  oscillations  were  stable  and  traveled  with  the  reflected  shock 
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Figure  3:  Comparison  of  pre-reflection  shock  wave  simulations  using  M 
(fine-dashed  line),  and  FCT  (coarse -dashed  line)  for  artificial 
diffusion.  Exact  solution  is  the  solid  line. 
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wave.  Monaghan-Gingold' s  viscosity  produced  no  oscillations,  however 
the  reflected  density  shock  front  was  more  diffused.  This  is  attributed 
to  the  reflected  shock  propagating  into  a  region  with  higher  density. 
There  are  more  downstream  particles  within  the  2h  smoothing  length  from 
the  shock  front  than  were  in  the  corresponding  region  for  the  pre- 
ref  lection  shock  wave.  The  imaginary  particle  method  compares  well  with 
the  exact  solution. 

Both  the  Lennard-Jones  and  spring  repulsive  force  schemes  produced 
essentially  indistinguishable  results.  Both  schemes  experienced 
particle  crossing  difficulties.  Figures  6  and  7  show  density,  velocity, 
temperature,  and  pressure  predictions  obtained  using  the  spring  force 
method  for  M-G  and  FCT  treatments,  respectively.  Flux  corrected 
transport  again  shows  a  tendency  to  produce  nonphysical  oscillations  on 
the  trailing  edge  of  the  reflected  shock  waves.  The  velocity  plots  of 
both  figures  show  instabilities  behind  the  reflected  shock  wave.  This 
is  due  to  particles  crossing  each  other  (though  not  the  end-wall)  near 
the  edge  of  the  applied  force.  The  instabilities  eventually  damp  out  as 
the  reflected  wave  propagates  away  from  the  tube  wall . 

The  density,  temperature,  and  pressure  plots  show  a  shortcoming  of 
the  repulsive  force  technique;  a  drop  in  value  along  the  boundary  edge. 
This  is  caused  by  the  particles  near  the  tube  end  missing  part  of  their 
summation.  The  spurious  pressure  and  density  drops  cause  an  initial 
velocity  toward  the  tube  end.  While  the  repulsive  force  acts  to  prevent 
this  initial  acceleration,  it  is  does  not  effect  the  values  of  the 
density  or  pressure,  and  they  remain  underpredicted.  An  attempt  was 
made  to  use  the  alternative  conservation  of  mass  equation  (Eq.  11),  as 
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Figure  6:  Plots  of  a)  velocity,  b)  pressure,  c)  density,  and  d)  temperature 

1.65  milliseconds  into  the  simulation,  and  after  reflection. 
Simulation  uses  Monaghan -Gingold  artificial  viscosity  (dashed 
line),  and  the  repulsive  force  (spring)  boundary  treatment. 
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this  method  calculates  the  density  using  the  velocity  divergence  instead 
of  particle  mass  summation.  Simulations  showed  that  after  the 
reflection  there  still  remained  a  drop  in  the  expected  value  because  of 
the  missing  parts  of  the  divergence  summation. 

3.3  Conclusion 

Conclusions  about  the  relative  performance  of  SPH  simulations  of 
reflected  shockwaves  using  artificial  diffusion  schemes  are  somewhat 
ambiguous.  Monaghan  and  Gingold's  viscosity  [9]  and  flux  corrected 
transport  [20]  have  been  compared.  Though  FCT  reproduces  the  sharpest 
shock  fronts,  it  appears  less  stable  after  the  reflection  and  more 
limited  in  application  than  M-G  viscosity.  However,  post-reflection 
simulations  using  M-G  viscosity  produce  more  heavily  smeared  shock 
waves. 

Conclusions  pertaining  to  different  boundary  treatments  are  less 
ambiguous.  The  most  accurate  simulations  of  the  reflected  shock  were 
made  using  the  imaginary  particle  scheme;  although  CPU  time,  memory 
requirements,  and  coding  complications  were  greater  than  for  other 
candidate  schemes.  Coding  techniques  can  reduce  the  extra  CPU  and 
memory  requirements  for  applications  in  higher  dimensions.  The 
repulsive  force  technique  is  the  easiest  to  code  and  requires  less  CPU 
time,  however  accuracy  is  sacrificed  and  particles  may  cross.  Neither 
the  immobile  particle  technique,  nor  Campbell's  [16]  pressure 
modification  boundary  technique  were  suitable  for  this  test  case. 
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Chapter  4:  INCOMPRESSIBLE  FLOW 

4-1  Introduction 

The  emphasis  of  the  previous  chapters  was  on  developing  the  SPH 
formalism  and  applying  the  method  to  one -dimensional ,  compressible, 
inviscid  flow.  The  object  of  this  chapter  is  to  extend  the  SPH  method 
to  model  two-dimensional,  viscous,  incompressible  flow.  The  test 
problem  used  to  accomplish  this  objective  are  the  shear  driven  cavity, 
Couette  flow,  and  Stoke 's  first  problem. 

4.2  Modifications 

The  standard  SPH  technique  for  solving  the  equations  of  motion 
utilize  summations,  based  on  particle  locations  and  a  weighting  kernel, 
to  determine  the  density,  velocity  and  internal  energy.  The  pressure  is 
a  thermodynamic  property  generally  obtained  from  the  density,  internal 
energy,  and  an  appropriate  equation  of  state. 

The  governing  equations  for  incompressible  viscous  flow  are: 


5  p/3 1  1  artificial  +  —  0 

(28) 

dpv/dt  =  -V*P  +  pH' 

(29) 

When  modeling  incompressible  flow,  the  standard  SPH  procedure  must 
be  modified  because  density  no  longer  effects  pressure  through  an 
equation  of  state.  This  modification  is  introduced  by  using  an 
artificial  compressibility  formulation. 
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4.2.1  Artificial  compressibility 

The  basis  for  the  artificial  compressibility  scheme  lies  in 
introducing  a  time  derivative  of  pressure.  This  term  can  be 
incorporated  into  an  equation  of  state .  The  following  equation  of  state 
for  water  has  been  applied  by  Monaghan  [14]  for  modeling  free  surface 
flows : 

P  -  B((p/po)^-l),  (30) 

with  r  =  7.415.  The  constant  B  is  simulation  specific,  and  depends  on 
the  maximum  fluid  velocity.  In  the  artificial  incompressibility 
technique,  the  speed  of  sound  is  defined  as  being  very  much  larger  than 
the  maximum  possible  speed  of  the  bulk  flow  (c  >  O.l^Vmax).  This  puts  a 
limit  on  the  maximum  change  in  the  density  values.  The  constant  B  is 
used  to  keep  the  density  changes  to  within  approximately  5%  of  the 
original  values.  In  these  simulations  B  was  defined  as  follows: 

B  =  (V„^0.1)^W7'.  (31) 

Where  7'  is  the  fluid's  specific  heat  ratio.  This  technique  models  an 
incompressible  fluid  as  slightly  compressible. 

4.2.2  Viscous  term 

A  common  requirement  for  all  the  previously  cited  SPH  formulations 
is  to  include  an  artificial  viscosity  component  to  the  dampen  the  high 
frequency  oscillations  inherent  to  Lagrangian  simulation  techniques. 
However,  in  simulating  viscous,  recirculating  flows  of  engineering 
interest,  true  viscosity  must  be  used. 

The  viscous  term  of  the  Eulerian  momentum  equation  is  written  as 
n'  -  1/p  V«(i/'Vv).  (31) 
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While  this  can  be  substituted  directly  into  the  SPH  formulation,  it  has 
been  shown  [22,23]  that  the  second  derivative  is  extremely  sensitive  to 
particle  disorder.  To  avoid  this  problem  in  the  heat  diffusion 
equation,  an  integral  approximation  was  used  to  model  the  second 
derivative  [22,23].  Because  of  the  direct  analogy  between  heat 
diffusion  and  viscous  momentum  diffusion,  this  technique  has  been 
applied  here  for  the  viscous  term.  Through  application  of  this 
technique,  the  viscosity  term  becomes: 

n'  =  -Smj^(2i/'(Vi-Vj)(rij-V3Wij))  /  (32) 

where  u'  is  the  kinematic  viscosity,  and  pij  is  the  averaged  density. 

4.3  Boundary  Treatment 

Using  the  results  from  the  previous  chapter  as  a  guide,  boundaries 
were  simulated  using  a  hybrid  scheme  involving  both  immobile  particles 
and  repulsive  forces.  Immobile  particles  were  placed  surrounding  the 
required  geometry  with  the  boundary  particles  exerting  the  repulsive 
force .  The  hybrid  technique  was  used  to  prevent  the  near  boundary 
density  dropoffs  intrinsic  to  the  repulsive  force  treatment  while 
maintaining  the  ease  and  minimum  CPU  time  requirements  characteristic  of 
the  repulsive  force  technique.  Three  layers  of  immobile  particles  were 
necessary  to  provide  the  correct  values  for  the  edge  particles.  The 
repulsive  force  was  modeled  using  the  Lennard- Jones  type  force  (Eqn.  26) 
and  was  exerted  solely  by  the  first  layer  of  immobile  particles. 
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4.4  The  Test  Problems 

The  first  test  problem  for  the  viscous,  incompressible  formulation 
was  Couette  flow.  Couette  flow  is  steady  viscous  flow  between  two 
infinitely  long  parallel  plates  with  the  lower  plate  stationary  and  the 
top  plate  moving  at  a  constant  speed.  The  fully  developed  solution 
yields  a  linear  velocity  profile.  This  test  problem  was  used  to 
determine  if  the  viscous  boundary  conditions  were  being  accurately 
modeled. 

Stoke 's  first  problem  was  the  next  test  case.  This  problem 
involves  the  transient  flow  development  for  a  fluid  near  a  suddenly 
moving  infinitely  long  plate.  Results  are  compared  to  the  exact 
solution.  This  problem  was  used  to  determine  if  diffusion  is  being 
properly  modeled  in  the  transient  case,  as  pressure  gradients  do  not 
contribute  significantly  to  the  fluid  motion  in  this  problem. 

The  third  test  problem  is  the  classic  shear  driven  cavity  flow. 
With  this  problem,  the  moving  top  of  a  closed  square  cavity  produces  a 
recirculation  pattern  in  the  fluid.  A  solution  was  obtained  for  a 
Reynolds  number  (UtopH/i/)  of  10.  Results  are  compared  with  solutions 
obtained  from  the  Eulerian  code  described  by  Miller  and  Schmidt  [24]. 
This  test  problem  tests  all  three  major  areas  of  interest:  the  boundary 
conditions,  the  viscous  diffusion,  and  the  artificial  compressibility. 

In  all  test  problems  the  transient  equations  were  stepped  in  time 
with  the  same  Runge-Kutta  scheme  described  previously.  No  artificial 
diffusion  has  been  included  in  any  of  the  simulations.  The  test 
problems  were  simulated  with  676  particles.  Of  these,  277  particles 
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were  used  for  the  boundaries,  and  were  not  independently  cycled  through 
the  summations . 

4.5  Discussion 

Figure  8  is  a  plot  of  the  velocity  profile  for  the  Couette  flow 
simulation.  This  result  matches  with  the  expected  solution  (constant 
slope).  The  velocity  profile  was  taken  at  a  time  well  after  the  steady 
state  solution  was  achieved.  This  indicates  that  the  hybrid  boundary 
condition  is  accurately  modeling  the  moving  surface,  and  the  viscous 
effects  are  being  properly  modeled  for  the  fully  developed  state. 

Figure  9  is  a  plot  of  the  transient  velocity  profile  for  Stoke 's 
first  problem.  The  SPH  solution  is  compared  with  a  similarity  solution 
[25].  Smoothed  particle  hydrodynamics  underpredicts  the  exact  solution 
by  approximately  10  percent  for  the  particle  closest  to  the  boundary. 
Particles  further  from  the  boundary  show  better  comparisons  with  the 
expected  solution.  It  is  expected  that  the  relatively  high  error  for 
the  particle  nearest  the  boundary  can  be  diminished  by  the  addition  of 
more  particles  near  the  moving  plate. 

Figure  10  shows  the  steady-state  streamlines  for  the  particles  in 
the  driven  cavity  flow  problem.  Figure  11  shows  the  actual  particle 
locations  for  the  same  time  (120  seconds).  These  plots  show  that  SPH 
reproduces  the  expected  physics  of  the  problem  and  the  particles 
maintain  an  essentially  equal  separation.  A  noticeable  error  lies  in 
the  upper  right  corner  of  the  cavity  in  figure  11.  An  area  depleted  of 
particles  form  soon  after  the  simulation  is  started.  Figures  12  and  13 
show  the  horizontal  and  vertical  velocity  profile  for  the  SPH 
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Figure  9:  Plot  of  transient  velocity  profile  for  Stoke 's  first  problem 
at  0.001  seconds  (dashed  line). 
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Figure  12:  Velocity  profile  along  vertical  centerline  for  the  driven 
cavity  with  a  Reynold's  number  of  10. 


36 


simulation,  respectively.  The  SPH  solutions  (dashed  lines),  are 
compared  with  the  Eularian  solutions  of  Miller  and  Schmidt  [24]  (solid 
lines).  Since  SPH  is  a  gridless  technique,  interpolation  was  necessary 
to  determine  the  velocity  profile.  Both  plots  show  an  underprediction 
in  the  velocity  values.  Because  the  Couette  flow  and  Stoke 's  first 
problem  simulations  compared  well  with  theory,  the  discrepancy  in  the 
driven  cavity  simulation  is  not  attributed  to  the  viscous  stress  or  the 
no -si ip  boundary  condition.  Rather,  the  error  is  attributed  to  the 
occurrence  of  the  particle  depletion  regions.  It  is  expected  that  with 
the  addition  of  more  particles  the  solution  will  converge  closer  to  the 
expected  result. 

4.6  Concliisions 

Conclusions  regarding  the  SPH  treatment  of  incompressible,  viscous 
flows  are  encouraging.  Smoothed  particle  hydrodynamics  successfully 
treats  the  Couette  flow  and  Stoke 's  first  problem.  This  indicates  that 
both  the  viscous  stress  and  the  no-slip  boundary  are  being  properly 
modeled.  The  driven  cavity  underpredicts  the  expected  solution  and  the 
error  is  attributed  to  a  loss  in  resolution  due  to  low  particle  areas. 


37 


Chapter  5 .  CONCLUSIONS/RECOMMENDATIONS 

Smoothed  particle  hydrodynamics  has  been  shown  to  have  the 
potential  to  accurately  model  flows  beyond  its  original  astrophysical 
applications.  Boundaries  have  been  implemented  successfully  through  a 
variety  of  techniques.  The  final  choice  of  boundary  technique, 
ultimately  lies  in  the  memory  and  CPU  requirements  of  the  user.  Flux 
corrected  transport  can  be  introduced  into  the  SPH  formulation  and  it 
has  advantages  and  disadvantages. 

Smoothed  particle  hydrodynamics  has  been  shown  to  have  great 
potential  in  modeling  incompressible  flow.  The  artifibial 
compressibility  treatment  and  the  introduction  of  true  viscosity 
reproduce  the  physics  of  the  test  problems,  though  not  to  a  desirable 
accuracy  in  the  driven  cavity.  The  area  of  error  has  been  determined  to 
be  a  loss  of  resolution  due  to  too  few  particles. 

Recommendations  for  further  study  include:  application  of  the 
artificial  diffusion  technique,  flux  corrected  transport,  to  higher 
dimensions;  finding  the  optimum  technique  for  introducing  viscosity; 
introducing  particles  in  depleted  regions ;  and  extending  the  technique 
to  modeling  multi-fluid  flows. 
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